##########################################################################################

library(data.table)
library(optparse)
library(WGCNA)
library(parallel)
library(dplyr)
library(RColorBrewer)

##########################################################################################

option_list <- list(
    make_option(c("--sample_list_file"), type = "character"),
    make_option(c("--rsem_file"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){

    sample_list_file <- "~/20220915_gastric_multiple/rna_combine/analysis/config/tumor_normal.list"
    rsem_file <- "~/20220915_gastric_multiple/rna_combine/analysis/images/wgcnv/CombineCounts.FilterLowExpression-MergeMutiSample.varianceStabilizingTransformation.tsv"
    out_path <- "~/20220915_gastric_multiple/rna_combine/analysis/images/wgcnv"

}

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

sample_list_file <- opt$sample_list_file
rsem_file <- opt$rsem_file
out_path <- opt$out_path

dir.create(out_path , recursive = T)

##########################################################################################
# The following setting is important, do not omit.
options(stringsAsFactors = FALSE);
# Allow multi-threading within WGCNA. This helps speed up certain calculations.
# At present this call is necessary.
# Any error here may be ignored but you may want to update WGCNA if you see one.
# Caution: skip this line if you run RStudio or other third-party R environments.
# See note above.
# enableWGCNAThreads()
# Load the data saved in the first part


##########################################################################################

info <- data.frame(fread(sample_list_file))
dat_tpm <- data.frame(fread(rsem_file))

##########################################################################################
## https://mp.weixin.qq.com/s/2VXprzgJSYJ39AhZ_FprCQ
## 1至少20个样本以上，越多越好；
## 2可以过滤点低表达或者低方差的基因，以减少干扰信息。但不太建议直接使用差异基因。
## 3WGCNA最初用于芯片测序数据，也适用于RNA-seq数据。关于RNAseq标准化，由于不涉及到不同基因之间的比较，所有常规标准化方式都可以。
## https://www.homedt.net/233619.html
## 输入数据形式如果有批次效应,需要先进行去除；
## 处理RNAseq数据，需要采用DESeq2的varianceStabilizingTransformation方法，或将基因标准化后的数据（如FPKM、CPM等）进行log2(x+1)转化
## https://www.biostars.org/p/280650/

##########################################################################################
##  选择合适的软阈值β
exp_dat <- dat_tpm[,1:(ncol(dat_tpm)-1)]
rownames(exp_dat) <- dat_tpm$gene_id
exp_dat <- t(exp_dat)

powers = c(c(1:10), seq(from = 12, to=20, by=2))
sft = pickSoftThreshold(exp_dat, powerVector = powers, verbose = 5)

out_name <- paste0( out_path , "/sft_plot.pdf" )

pdf(out_name)
cex1 = 0.9
par(mfrow = c(2,1))
### （1）是否符合幂律分布；
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
     main = paste("Scale independence"));
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     labels=powers,cex=cex1,col="red")
abline(h=0.80,col="red")
### （2）节点的平均连接度
plot(sft$fitIndices[,1], sft$fitIndices[,5],
     xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
     main = paste("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")
#par(mfrow = c(1,1))
dev.off()

## 基本在6时，拟合幂律分布的结果是比较好的；同时节点的凭据连通性也趋于稳定了。之后会使用这个参数建立网络
## To maximize strong correlations between genes, we assigned a soft-thresholding power β = 8 since it was the 
## lowest β with the highest R2 value that passed the 0.8 threshold 
power <- sft$fitIndices[,1][-sign(sft$fitIndices[,3])*sft$fitIndices[,2] > 0.8][1]
#power <- 6

##########################################################################################
## 建立网络，鉴定模块
## https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/Tutorials/FemaleLiver-02-networkConstr-blockwise.pdf
tom_file <- paste0( out_path , "/TOM-blockwise" )

## I am not aware of a principle from which one could derive an “appropriate” value, but in practice, 
## on data sets with 50-100 samples, using 0.15 to 0.2 has worked fairly well. 
## For fewer samples a larger valuse (0.25 to 0.3) may be warranted. 
## If you want larger modules, increase the value; 
## if you want smaller modules at the risk of having redundant modules 
## (modules with very similar functional annotation and very similar fuzzy module membership), 
## you can decrease the value to say 0.10, maybe even lower if you have lots of samples (hundreds or more).

set.seed(42)
net = blockwiseModules(exp_dat, 
     power = power , TOMType = "unsigned", minModuleSize = 100,
     numericLabels = TRUE)

MEs = net$MEs # 模块特征向量
#TOM矩阵层次聚类结果
geneTree = net$dendrograms[[1]]

out_file <- paste0( out_path , "/plotDendroAndColors.pdf" )
sizeGrWindow(12, 9)
pdf(out_file)
mergedColors = labels2colors(net$colors) 
#pdf(file = 'D:/practice/test/module.pdf',width = 12,height = 8)
plotDendroAndColors(net$dendrograms[[1]],
                  mergedColors[net$blockGenes[[1]]], 
                  "Module colors",
                  dendroLabels = FALSE,
                  hang = 0.03,
                  addGuide = TRUE,
                  guideHang = 0.05)
dev.off()



moduleColors <- labels2colors(net$colors)
# Recalculate MEs with color labels
MEs0 = moduleEigengenes(exp_dat, moduleColors)$eigengenes
# 计算根据模块特征向量基因计算模块相异度：
MEDiss = 1 - cor(MEs0);
# Cluster module eigengenes
METree = hclust(as.dist(MEDiss), method = "average");

out_name <- paste0( out_path , "/Clustering_module_eigengenes.pdf" )
pdf(out_name)
plot(METree,
     main = "Clustering of module eigengenes",
     xlab = "",
     sub = "")
# 在聚类图中画出剪切线
MEDissThres = 0.25
abline(h = MEDissThres, col = "red")
dev.off()


## power：软阈值的选择
## corType：计算相关性的方法；可选pearson(默认)，bicor。后者更能考虑离群点的影响。
## networkType:计算邻接矩阵时，是否考虑正负相关性；默认为"unsigned",可选"signed", "signed hybrid"
## TOMType：计算TOM矩阵时，是否考虑正负相关性；默认为"signed",可选"unsigned"。但是根据幂律转换的邻接矩阵(权重)的非负性，所以认为这里选择"signed"也没有太多的意义。
## minModuleSize：模块的最少基因数
## mergeCutHeight：合并模块的阈值
## numericLabels：模块名是否为数字；若设置FALSE，表示映射为颜色名。
## saveTOMs：是否保存TOM矩阵；如果设为TRUE，需要设置saveTOMFileBase参数，提供保存文件名；设置numericLabels参数，是否将模块名保存为颜色名
## nThreads：交代线程数，适用于Linux环境
## verbose：0默认安静的执行，值越大表示给出的运行提示信息越多。


##########################################################################################
## 2.5 挑选模块Hub基因
## (1) 模块内基因连接度
# WGCNA offers a function "intramodularConnectivity" that actually also computes the total connectivity (within the whole network). 
# In addition to total and intramodular connectivity it also gives the extra-modular connectivity (kTotal-kWithin), 
# and the difference of the intra- and extra-modular connectivities (kIn-kOut) for all genes. 
# The value returned from the function is a data frame with one column for each of the measures mentioned above. 
# Intramodular connectivity is calculated by summing adjacency entries (excluding the diagonal) to other nodes within the same module. 
# Intramodular connectivity (kWithin)
# After raising the module membership to a power of 6, it is highly correlated with the intramodular connectivity (kWithin).

adjacency = adjacency(exp_dat, power = power)
# TOM = TOMsimilarity(adjacency)
Alldegrees <- intramodularConnectivity(adjacency, net$colors)

dat_cluster <- data.frame( cluster = as.numeric(net$colors) )
rownames(dat_cluster) <- names(net$colors)
dat_cluster <- cbind( Alldegrees , dat_cluster )

#### kTotal:基因在整个网络中的连接度
#### kWithin: 基因在所属模块中的连接度，即Intramodular connectivity
#### kOut: kTotal-kWithin
#### kDiff: kIn-kOut
out_name <- paste0(out_path , "/wgcnv_HubModule.tsv")

## 模块所在连接度前5%的基因作为hub基因
## https://jneuroinflammation.biomedcentral.com/articles/10.1186/s12974-019-1433-4
result <- c()
## 0代表没有聚成一类
for(cls in unique(dat_cluster$cluster)[unique(dat_cluster$cluster)!=0]){
    tmp <- subset( dat_cluster , cluster == cls)
    t_th <- as.numeric(quantile(tmp$kWithin , seq(0,1,0.05))[20])

    tmp$Hub <- ifelse(tmp$kWithin >= t_th , "TRUE" , "FALSE")
    tmp <- tmp[order(tmp$kWithin , decreasing=T),]
    tmp$kWithin_rank <- order(tmp$kWithin , decreasing=T)/nrow(tmp)
    result <- rbind(result , tmp)
}

result$gene_id <- rownames(result)

write.table( result , out_name , row.names = F , sep = "\t" , quote = F )

#### https://github.com/halryd/high_S100A_hub_genes/blob/master/08_chap3_select_hub_genes.Rmd
## In the tutorials of WGCNA they indicate three different way to measure the degree to which a gene is a hub. 
## It is shownd that GS and MM correlate quit well an that. MM and 
## 1. Gene significance (GS)
## 2. Module membership (MM)
## 3. Module connectivity (MC?)


## - In paper "Inflammatory, regulatory, and autophagy co-expression modules and hub genes underlie the peripheral immune response 
## to human intracerebral hemorrhage" they are using top 5% highest membership to their respective module as definition of a hub gene-
## References for 5% membership: @RN7


## In "Gene co-expression network analysis reveals common system-level properties of prognostic genes across cancer types" they  defined the 1% (or 5%) of nodes with the highest connectivity as hubs 25–27.
## They write further "Given a network, we then obtained several key network properties such as the edge weight, node connectivity and modularity. Connectivity was defined as the sum of the weights across all the edges of a node, and the top 1% (or 5%) of the genes with the highest connectivity in the network were defined as hub genes."
## References for 5% connectivity @RN5

